; docformat = 'rst'
;
; NAME:
;   MR_LITTLEREGION
;
;+
; :Purpose:
;   Find values in array surrounding a particular grid cell
;
; :Categories:
;   Model utilities
;
; :Returns:
;   Array of values surrounding a grid cell.  If regionsize eq 4 or 8 returns a 1D array (see below).
;	  If RegionSize>9 returns square array of side sqrt(regionsize).
;
; :Params:
;   Array : 2D array
;		LonI : Array index in first dimension
;		LatI : Array index in second dimension
;   
; :Keywords:
;   RegionSize : Integer specifiying number of surrounding grid cells to return.
;			Allowed values::
;				4 (DEFAULT) returns cell above, right, below, left, in 1D array in that order
;				8: returns 1D array of 8 grid cells surrounding, starting above and then clockwise
;				9, 25, 49... square of odd-numbered side length: Returns 2D square array of values
;					centered on LonI, LatI.
;   Indices : (output) the column and row indices in Array corresponding to the output values
;   NoWrap : (input) do not wrap around.  Default is to wrap in first dimension (e.g. longitude)
;
; :Examples:
;   Find 3x3 region around grid cell at position 2, 2, in matrix a::
;		   IDL> a=findgen(5, 5)                
;		   IDL> print, a                                             
;         0.00000      1.00000      2.00000      3.00000      4.00000
;         5.00000      6.00000      7.00000      8.00000      9.00000
;         10.0000      11.0000      12.0000      13.0000      14.0000
;         15.0000      16.0000      17.0000      18.0000      19.0000
;         20.0000      21.0000      22.0000      23.0000      24.0000
;		   IDL> print, mr_littleregion(a, 2, 2, regionsize=9, indices=i)
;         6.00000      7.00000      8.00000
;         11.0000      12.0000      13.0000
;         16.0000      17.0000      18.0000
;		   IDL> print, i                                             
;           1           1
;           2           1
;           3           1
;           1           2
;           2           2
;           3           2
;           1           3
;           2           3
;           3           3
;
; :History:
;   Written by: Matt Rigby, MIT, May 31st 2011
;
;-
Function mr_littleregion, array, lonI, latI, regionsize=regionsize, indices=indices, nowrap=nowrap

  compile_opt idl2

  lonsize=n_elements(array[*, 0])
  latsize=n_elements(array[0, *])

  if keyword_set(regionsize) eq 0 then regionsize=4

  if LatI eq 0 then latIM=latI else latIM=latI-1
  if latI eq latSize-1 then latIP=latI else latIP=latI+1
  
  if lonI eq 0 then lonIM=LonSize-1 else lonIM=lonI-1
  if lonI eq LonSize-1 then lonIP=0 else lonIP=lonI+1

  if keyword_set(nowrap) then begin
    if lonI eq 0 then lonIM=0 else lonIM=lonI-1
    if lonI eq LonSize-1 then lonIP=LonSize-1 else lonIP=lonI+1
  endif
  
  case regionsize of
  4: begin
    littleregion=[array[lonI, latIP], $
      array[lonIP, latI], $
      array[lonI, latIm], $
      array[lonIM, latI]]
    indices=[ [lonI, latIP], [lonIP, latI], [lonI, latIM], [lonIM, latI] ]
    end
  8: begin
    littleregion=[array[lonI, latIP], $
      array[lonIP, latIP], $
      array[lonIP, latI], $
      array[lonIP, latIM], $
      array[lonI, latIM], $
      array[lonIM, latIM], $
      array[lonIM, latI], $
      array[lonIM, latIP]]
    indices=[[lonI, latIP], $
      [lonIP, latIP], $
      [lonIP, latI], $
      [lonIP, latIM], $
      [lonI, latIM], $
      [lonIM, latIM], $
      [lonIM, latI], $
      [lonIM, latIP]]
    end

  else: begin
    side=round(sqrt(regionsize), /l64)
    latM=max([0, latI - side/2])
    latP=min([latsize-1, latI + side/2])
    littleregion=dblarr(side, LatP - LatM+1)

    indices_temp=indgen(LonSize, LatSize, /long)
    
    if keyword_set(nowrap) then begin
      lonM=max([0, lonI - side/2])
      lonP=min([lonsize-1, lonI + side/2])
      littleregion=array[LonM:LonP, LatM:LatP]

      indices=indices_temp[LonM:LonP, LatM:LatP]
    endif else begin
      if lonI ge side/2 and lonI le lonsize-side/2-1 then begin
        littleregion=array[lonI - side/2:lonI + side/2, latM:latP]
        indices=indices_temp[lonI - side/2:lonI + side/2, latM:latP]
      endif else begin
        indices=lonarr(side, LatP - LatM+1)
        if lonI lt side/2 then begin
          littleregion[0:side/2 - LonI - 1, *]=array[LonSize - side/2 + LonI:LonSize-1, LatM:latP]
          littleregion[side/2 - LonI:side-1, *]=array[0:side/2 + lonI, LatM:LatP]
          indices[0:side/2 - LonI - 1, *]=indices_temp[LonSize - side/2 + LonI:LonSize-1, LatM:latP]
          indices[side/2 - LonI:side-1, *]=indices_temp[0:side/2 + lonI, LatM:LatP]
        endif
        if lonI gt lonSize - side/2-1 then begin
          littleregion[0:side/2 + lonsize-1 - lonI, *]=array[LonSize - (lonsize-lonI) - side/2:LonSIze-1, LatM:LatP]
          littleregion[side/2 + lonsize - lonI :side-1, *]=array[0:side/2 - (lonsize - lonI), latM:LatP]
          indices[0:side/2 + lonsize-1 - lonI, *]=indices_temp[LonSize - (lonsize-lonI) - side/2:LonSIze-1, LatM:LatP]
          indices[side/2 + lonsize - lonI :side-1, *]=indices_temp[0:side/2 - (lonsize - lonI), latM:LatP]
        endif
      Endelse
      
    endelse

    ncol=n_elements(indices[*, 0])
    nrow=n_elements(indices[0, *])
    ai=lonarr(n_elements(indices))
    for n=0, nrow-1 do ai[n*ncol:(n+1)*ncol-1]=indices[*, n]

    indices=array_indices(array, ai)
   end
  EndCase

  return, littleregion

End
